Coulomb interaction and transient charging of excited states in open nanosystems 
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We obtain and analyze the effect of electron-electron Coulomb interaction on the time dependent 
current flowing through a mesoscopic system connected to biased semi-infinite leads. We assume 
the contact is gradually switched on in time and we calculate the time dependent reduced density 
operator of the sample using the generalized master equation. The many-electron states (MES) of 
the isolated sample are derived with the exact diagonalization method. The chemical potentials of 
the two leads create a bias window which determines which MES are relevant to the charging and 
discharging of the sample and to the currents, during the transient or steady states. We discuss the 
contribution of the MES with fixed number of electrons N and we find that in the transient regime 
there are excited states more active than the ground state even for N — 1. This is a dynamical 
signature of the Coulomb blockade phenomenon. We discuss numerical results for three sample 
models: short ID chain, 2D lattice, and 2D parabolic quantum wire. 

PACS numbers: 73.23.Hk, 85.35.Ds, 85.35.Be, 73.21. La 



I. INTRODUCTION 

Due to the increasing interest in ultra-fast electron 
dynamics considerable progress occurred recently in the 
theoretical description of time dependent mesoscopic 
transport. New methods and numerical implementations 
are rapidly evolving. Transient currents in open n anos- 
tructures are studied with Green-Keldysh formalism^^ 
scattering theory,"^ and quantum master equation.^^^^^^ 
Most of the results were obtained for noninteracting elec- 
trons due to the well known computational difficulties to 
include time-dependent Coulomb effects. 

It is nevertheless clear that the electron-electron inter- 
action is important in such problems. An effort to incor- 
porate it has been recently done by Kurth et al^ followed 
by Myohanen et a/J^who have described correlated time- 
dependent transport in a short ID chain defined by a 
lattice Hamiltonian. The ID sample was connected to 
external leads and the current was driven by a time- 
dependent bias. Those authors used a method based 
on the Kadanoff-Baym equation for the non-equilibrium 
Green's function combined with the time-dependent den- 
sity functional theory to include the Coulomb interac- 
tion in the sample. Once the Green's functions were 
calculated total average quantities of interest could be 
obtained, like charge density or current, both in the tran- 
sitory and in the steady state. However this method does 
not say much about the dynamics of specific internal 
states of the saniple system. In view of the spectroscopy 
of excited statea^I' it is important to have a theoretical 
tool for understanding separately the charging and re- 
laxation of the ground states and excited states in meso- 
scopic systems in time-dependent conditions. 

Our alternative is to use the statistical, or density op- 
erator. The complete information about the time evo- 
lution of each quantum state of the sample is captured 
in the reduced density operator (RDO), which is the so- 



lution of the generalized master equation (GME). Once 
the RDO is defined in the Fock space the inclusion of 
the Coulomb interaction becomes a known computational 
problem: obtaining the many-electron states (MES) of 
the sample. The RDO matrix is then calculated in the 
basis of the interacting MES. 

Let us enumerate some of the previous theoretical 
schemes to treat transport and electron-electron inter- 
action with the master equation. One of the first at- 
tempts to derive a master equation for an interacting sys- 
tem with time-dependent perturbations belongs to Lan- 
greth and Nordlander for the Anderson model. Gurvitz 
and Prager started from the time-dependent Shrodinger 
equation for the MES wave functions and ended up with 
Bloch-like rate equations for the density matrix of a quan- 
tum dot.^^ The electronic currents were calculated in the 
steady state and it was shown that the Coulomb interac- 
tion renormalizes the tunneling rates between the leads 
and the system. In the same context Konig et a/.^^ de- 
veloped a powerful diagrammatic technique by expand- 
ing the RDO of a mesoscopic system in powers of the 
tunneling Hamiltonian. The time-dependence of the sta- 
tistical operator of the coupled and interacting system 
implies a quantum master equation for the so called pop- 
ulations. In this method the Coulomb interactions are 
treated exactly, which makes it appealing for studying 
various correlation effects like cotunneling.tl^ The con- 
nection between the real-time diagrammatic approac h of 
Konig et al}^ and the Nakajima-Zwanzig approach 
to the generalized master equation (GME) approach was 
made transparent by Timm.^^ 

More recently Li and Yan^^ combined the n-resolved 
master equation and the time dependent density- 
functional method to write down a Kohn-Sham master 
equation for the reduced single-particle density matrix. 
Also, Esposito and Galperin,^ using the equation of mo- 
tion for the Hubbard operators, have obtained a many- 
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body description of quantum transport in an open sys- 
tem and established a connection between the GME and 
non-equihbrium Greeen's functions. They studied simple 
systems in the steady state regime: a resonant level cou- 
pled to a a single vibration mode, an interacting dot with 
two spins, and a two-level bridge. Another recent work 
by Darau et al?^ implemented the GME for a benzene 
single-electron transistor and used exact MES to compute 
steady state currents within the Markov approximation.^ 
The stability diagram and the conductance peaks were 
obtained and a current blocking due to interferences be- 
tween degenerated orbitals was noticed. 

In our previous paper^^ we considered the GME 
method for the RDO of independent electrons in the Fock 
space. We discussed the transient transport through 
quantum dots and quantum wires. The contact between 
the leads and the sample was switched on at a certain ini- 
tial moment t^. We discussed extensively the occupation 
of the states within the bias window and the geometrical 
effects on the transient currents. We described the cou- 
pling between the sample and the leads via a tunneling 
Hamiltonian in which we took into account the spatial 
extension of the wave functions of both subsystems in 
the contact region. 

In spite of earlier or more recent attempts a complete 
description of the Coulomb effects in the time-dependent 
transport is still missing, especially in sample models 
larger than a few sites. In the present work we com- 
bine the GME method with the Coulomb interaction in 
the sample and we analyze the dynamics of the electrons 
starting with the moment when the leads are coupled to 
the sample until a steady state is reached. The Coulomb 
interaction is included in the Hamiltonian of the isolated 
sample and the interacting MES are calculated with the 
exact diagonalization method. This means the Coulomb 
interaction is fully included with no mean field assump- 
tion or density- functional model. The number of single- 
electron states (SES) used to define the matrix elements 
of the Hamiltonian of interacting electrons is sufficiently 
large such that the MES of interest are convergent. Due 
to the finite bias window only a limited number of MES 
participate to the charge transport through the sample, 
i. e. only those energetically compatible with the elec- 
trons in the leads. Hence the MES of interest are selected 
by the chemical potentials in the leads. We calculate the 
RDO matrix elements in the subspace of these MES using 
the GME. The electron-electron interaction in the leads 
is neglected. 

It is well known that the Fock space increases expo- 
nentially with the number of SES. In addition the time 
dependent numerical solution of the GME is also com- 
putational expensive. So at this stage we are limited to 
describe only few electrons in the system: up to five in a 
small system, but only up to three in a larger one. 

The paper is organized as follows. In Section 2 we 
briefly describe the GME, the inclusion of the Coulomb 
interaction, and the selection of the MES. Next, in Sec- 
tion 3, we show results for three models: a short ID 



chain, a 2D lattice of 12x10 sites, and a finite quantum 
wire with parabolic lateral confinement. Conclusions and 
discussions are presented in Section 4. 

II. GME METHOD AND COULOMB 
INTERACTION 

In this section we summarize the main lines of our 
method. The equations apply both to the lattice and 
continuous models. The time-dependent transport prob- 
lem is considered within the partitioning approach which 
is known both from the pioneering work of Caroli^^ and 
from the derivation of the GME. Prior to an initial time 
to the left lead (L) having a "source" role, and the right 
lead (R) having a "drain" role, are not connected to the 
sample and therefore can be characterized by equilibrium 
states with chemical potentials fiL and jiR respectively. 
Our aim is to compute the time dependent currents flow- 
ing through the sample and leads starting at moment to, 
when the three subsystems are connected, until a station- 
ary state is reached. 

The generic Hamiltonian of the total system consisting 
of the sample plus the leads is: 

H{t)=HL + Hn + Hs+HT{t). (1) 

Hi with I = L^R are the Hamiltonians of the leads. We 
denote by Sqi and ipqi the single-particle energies and 
wave functions respectively, for each lead. Using the cre- 
ation and annihilation operators associated to the single- 
particle states, cj^ and c^/, we can write 

Hi = J dqeqicl^Cqi . (2) 

Hs is the Hamiltonian of the sample. In the absence 
of the interaction the SES have discrete energies denoted 
as En and corresponding one-body wave functions (j)^ (r) . 
Using now the creation and annihilation operators for the 
sample SES, dj^ and d^, we can write 

Hs ^Endldn + ^ X] Vnm,n^m^dldl^dm^dn^ . (3) 
fi nm 

The second term in Eq. (|3| is the Coulomb interaction. 
In the SES basis the two-body matrix elements are given 
by: 

Vnm,n'm' = j drdr' (f>l{r)(p^{r')u{r - r' )</>»' W 0m' (r')> 

(4) 

where u{r — r') is the Coulomb potential. 

The third term of Eq. ([T]) is the so-called tunneling 
Hamiltonian describing the transfer of particles between 
the leads and the sample: 

Hrit) = E E / dqXiit)iTl„clid„ + h.c.) . (5) 

l=L,R n 
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Ht contains two important elements: (1) The time de- 
pendent switching functions xi if) which open the contact 
between the leads and the sample; these functions mimic 
the presence of a time dependent potential barrier. (2) 
The coupling T^^ between a state with momentum q of 
the lead / and the state n of the isolated sample, with 
wave function (/)^. The coupling coefficients T^^ depend 
on the energies of the coupled states and, maybe more 
important, on the amplitude of the wave functions in the 
contact region. As we have shown in our previous worli^ 
this construction allows us to capture geometrical effects 
in the electronic transfer. A precise definition of the cou- 
pling coefficients is however model specific, and will be 
mentioned in the next section. 

The evolution of our system is completely determined 
by the statistical operator W{t) associated to the total 
Hamiltonian H{t) defined in Eq.Q. W{t) is the solution 
of the quantum Liouville equation with a known initial 
value, prior to the coupling of the sample and leads: 

ihW{t) = [H{t), W{t)] , W{t < to) = PLPRPs , (6) 

The isolated leads are described by equilibrium distribu- 
tions, 

" Tr^{e-/5(^^-/^^^0} ' ^ = ^'^' (^) 

and the isolated sample by the density operator ps- Af- 
ter the coupling moment the dynamics of the sample is 
conveniently described by the RDO which is defined by 
averaging the total statistical operator over those degrees 
of freedom belonging to the leads: 

pit) = TiLTinWit), p{to) = PS- (8) 

In the absence of the electron-electron interaction the 
MES eigenvectors of Hs are bit-strings of the form \u) = 
\ii,i2^ ..,^5^...), where = 0,1 is the occupation number 
of the n-th SES. The set {u} is a basis in the Fock space 
of the isolated sample and the RDO can be seen as a 
matrix in this basis. From Eqs. (|6|-([8| we obtain in the 
lowest (2-nd) order in the coupling parameters T^^ the 
GME (see Ref. 7 for details): 

m= - {[Hs.Pd)] 

- ^ E j dqxim%u%m^h.c.),{^) 

where the coupling operator Tqi has matrix elements 

(^.0^. = E^«nH^V)- (10) 

n 

The operators Qqi and Uqi are defined as 

Qqiit) = e-^^^- f dsxi{s)Uqi{s)e'^'-'^'^^e''''', 

J to 



and fi is the Fermi function of the lead 

In the presence of the electron-electron interaction in 
the sample the MES which are eigenstates of Hs are lin- 
ear combinations of bit-strings: Hs\ol) = £a\o^), where 
I a) = ^^j^CauW)^ Cau being the mixing coefficients 
which can be found together with the energies Ea by 
diagonalizing Hs. (To distinguish better between the 
noninteracting and the interacting MES we use the right 
angular bracket for the former and the regular curved 
bracket for the later.) Using now the set {a} as a basis, 
i. e. the interacting MES, the GME has the same form 
as Eq. ([9|, where the matrix elements of all operators 
are now defined in the interacting basis and the matrix 
elements of the coupling operators are 

n 

Because the sample is open the number of electrons N 
contained in the sample is not fixed. The Hamiltonian 
Hs given in Eq. ([3| commutes with the total "number" 
operator dl^^dn. Thus is a "good quantum number" 
such that any state \a) has a fixed number of electrons. 
So the MES can also be labeled as \a) = \N,i) with 
i = 0,1,2,... an index for the ground and excited states 
of the MES subset with N electrons. The many-body 
energies can also be written as Sa = . In the practical 
calculations N varies between (the vacuum state) and 
J^max which is the total number of SES considered in the 
numerical diagonalization of Hs- The total number of 
MES is thus 2^— . 

If the coupling between the leads and the sample is 
not too strong we expect that only a limited number of 
MES participate effectively to the electronic transport. 
These states are naturally selected by the bias window 
[Pr^ Pl]- In the following examples, by selecting suitable 
values of the chemical potentials in the leads, we will 
truncate the basis of interacting MES to a reasonably 
small subset such that we can solve numerically Eq. ([9| 
with our available computing resources. To relate the 
bias window with the effective MES we need to consider 
the chemical potential of the isolated sample containing 
N electrons, 

4^ = 4) -4°!,, (12) 

which is the energy required to add the N-th electron on 
top of the ground state with A^ — 1 to obtain the i-th MES 
with A^ particles. We expect the current associated to 
the MES I A^, i) to depend on the location of the chemical 

potential p^^ relatively to the bias window. In particular 
it is clear that if at the coupling moment to the sample is 

(i) 

empty all MES with p}/ Pl will remain empty both 
during the transient and the steady states, so they can 
be safely ignored when solving the GME. 
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III. 



MODELS AND RESULTS 



We have numerically implemented the GME method 
both for lattice and continuous models. The sample mod- 
els are: a short ID chain with 5 sites, a 2D rectangular 
lattice with 12 x 10 = 120 sites, and a short quantum 
wire withe parabolic lateral confinement. In all cases the 
coupling functions have the form 



xiit) = i 



eT* + 1 



(13) 



with 7 a constant parameter, such that at the initial mo- 
ment, which is to = 0, we have Xi{^) = (no coupling), 
and in the steady state, for t ^ oo, = 1 (full coupling). 



A. A toy model: short ID chain 

In this model the two semi-infinite leads are attached 
to the ends of a ID chain with 5 sites. The coupling 
between a lead state with wave function ipqi and a sam- 
ple state with wave function (j)n is given by the product 
between the wave functions at the contact site: 



qn 



Virgimnin), 



where is the contact site of the lead / 
sites of the sample being il = 1 and ir - 



(14) 

= L, the end 
5. 
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FIG. 1: (Color online) The equilibrium chemical potentials 
/i^'* for 1 < A/" < 5 as a function of the interaction strength 
U . The dotted lines mark the chemical potentials of the leads 
selected in the transport simulations shown in the next figure, 
i.e. iiL — 5.25 and fiR — 4.75. 

The reason to call this a toy model is that we can ob- 
tain the complete set of 2^ = 32 MES, i. e. we do not 
need to cut the basis of the 5 SES. We also do not need 
to cut the MES basis, all matrix elements of the statisti- 
cal operator can be numerically calculated, even if not all 
of them might be important for the currents. In addition 
we will consider the strength of the Coulomb interaction 
as a free parameter /7, whereas in a realistic systems this 
is fixed by the electron charge and the dielectric con- 
stant of the material. Our goal is to have a qualitative 



understanding of the underlying physics, and in particu- 
lar to show the presence of the Coulomb blocking effects 
at certain values of U or of the chemical potentials of the 
leads. The Coulomb matrix elements defined in Eq. Q 
are calculated as 



(t>n'{i)(t>ni'{i') ' (15) 



h - ^- I 

In Fig.l we show the equilibrium chemical potentials 
/i^^ corresponding to ground states with 1 < A/" < 5 par- 
ticles against the interaction strength \J . One observes 
a linear dependence of /i^^ on [/, with slope increasing 
with N . Obviously the total Coulomb energy increases 
both with U and N. 

Let us now briefiy review the Coulomb blockade 
scenario Suppose the isolated sample contains N elec- 
trons and the chemical potentials of the leads are cho- 
sen such that < /iR < /iL < Mtv+i- Then the bias 
window [fiR^ jul] may include one or more of the excited 
configurations with N particles. In general some states 
with N electrons may have excitation energies exceed- 
ing /iL or even /i^^i- This situation corresponds to the 
Coulomb blockade phenomenon. Indeed, the addition of 
the {N + l)-th electron is energetically forbidden. Con- 
sequently the current in the steady state should vanish. 
However, shorter or longer transient currents are gener- 
ated by all many-body configurations in the vicinity of 
the bias window. 

Fig. 2(a) and 2(b) show the total currents in the left 
lead and the total charge residing in the sample for sev- 
eral values of the interaction strength. U is measured 
in units of ts, the hopping parameter in the sample,'^ 
and the time is expressed in units of h/ts while the cur- 
rent is in units of ets/h. The coupling constant in Eq. 
(13) is 7 = 1. The system is initially empty and thus 
p(0) = |00000)(00000|. 

The chemical potentials of the leads, fiL = 5.25 and 
jiiR = 4.75, are chosen such that in the absence of 
Coulomb interaction, i. e. for /7 = 0, /i4^^ is located 
within the bias window. In this case we obtain in the 
steady state the mean number of electrons about 3.6 and 
a non- vanishing current in the leads. This is understand- 
able, since /i^^"^ =£^4 = 5, which is the 4-th level of 
the isolated sample. The occupation of this level in the 
steady state is about 0.6, the other states being either full 
or empty. Also in this case, the excited states have small 
contributions to the steady state current as the system 
tends to be in the ground state with N = 3 electrons. 
Those contributions may also depend on the coupling 
strength of individual states with the leads, but in gen- 
eral remain small. 

The situation may change for U 7^ 0. For the inter- 
acting system, e. g. for U = 0.3, the system settles down 
in the Coulomb blockade regime, the total current be- 
ing almost suppressed in the steady state. This happens 
because the interaction pushes the chemical potentials 
upwards such that for U = 0.3 both ground states with 
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FIG. 2: (Color online) The total current entering the 5x1 
sample from the left lead as a function of time for the different 
values of the interaction strength U . The chemical potentials 
of the leads /xl = 5.25 and iir — 4.75. 



A/" = 3 and N = A electrons are outside the bias win- 
dow and cannot produce steady currents. When the in- 
teraction strength is further increased to U = 0.5 and 
U = 1 the steady state currents are gradually restored. 
This could look surprising, but one can see in Fig.l that 
by increasing U the ground state configuration with 3 
electrons approaches and enters the bias window. Con- 
sequently the transport becomes again possible. Note 
that while the steady state currents are not monotonous 
w.r.t. U the charge absorbed in the system continuously 
decreases, Fig. 2(b). 

In transport experiments the strength of the electron- 
electron interaction is indeed fixed. The usual way to 
obtain the Coulomb blockade is to vary the chemical po- 
tentials of the leads relatively to the energy levels of the 
sample, or vice versa. In Fig. 3 we show the currents 
in both leads for different values of the chemical poten- 
tial /ii^, while keeping fixed /iL = 7. The strength of 
the Coulomb interaction is /7 = 1 and fi^^^ almost equals 
/j^L. The steady state value of the current decreases as 
jj^R increases, because fewer states are included in the 
bias window. The Coulomb blockade onset occurs for 
fiR > 5, when /ig^^ drops below fiR. We observe that the 
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FIG. 3: The time-dependent total currents in the left and 
right leads at different values of the chemical potential /jr. 
The current in the right lead starts at negative values. Other 
parameters: Vl — Vr — 0.750, U = 1.0. 



maximum value of the total current in the left lead does 
not change much when /ir varies. In contrast, the tran- 
sient current in the right lead is negative and increases 
in magnitude as jur increases. This means that the right 
lead feeds the many-body configurations that fall below 
Mi?. 

The contribution of the excited states to the transient 
and steady state currents depends strongly on the bias 
window. In Fig. 4 we show the currents entering the sam- 
ple from the left lead, carried by the states with N = 2 
and = 3 electrons, for /i^^ = 3, 4, 5 (the cases with non- 
vanishing current in the steady state). We also show sep- 
arately the contribution to the currents associated to the 
ground state configurations, related to fi^"^ and fi^\ and 
the complementary contribution of all the excited states 
with 2 and 3 particles. In this case the wave vectors of 
the ground states are mostly given by the non-interacting 
wave vectors: 1 11000) with weight 97% and 1 11100) with 
98% for TV = 2 and TV = 3 respectively. 

For jj^R = 3 the steady state current of the ground 
state configuration is vanishingly small and so the total 
negative current associated to two-particle states comes 
mostly from the excited states. In the many-body energy 
spectrum of the isolated sample we obtain 5 excited con- 
figurations with /i2*^ G [/ii?,/iL] = [3,7]. As /J.R moves up 
the steady state current of the ground state with N = 2 
becomes also negative. The combined contributions of 
the excited states vanishes at /i^^ = 5. As can be seen 
from Fig. 1 /j^r = b is well above iJ.'^\ but very close to 

/jj^^ . Consequently, the ground configuration with N = 2 
is heavily populated in the steady state, whereas the ex- 
cited states have low probability and thus weak current. 
Actually, as we have checked, all the currents associated 
to each excited state with N = 2 vanish individually. In 
the transient regime however the N = 2 currents in all 
three cases are dominated by the excites states. 

The currents of the excited states having N = 3 elec- 
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FIG. 4: The separate contributions to the current of the 
ground state with N particles and of all excited states with 
N particles, for different values of /hr. For completeness we 
also include the total currents Jl for the same configura- 
tions. The discussion is made in the text. Other parameters: 
Vl = Vr = 0.750, U = 1.0. 



trons are positive at /ii? = 3, but change sign at /i^^ = 4. 
For fiji = 5 their magnitude exceeds the contribution of 
the ground state which is always positive. A more de- 
tailed analysis of the currents carried by specific excited 
states will be given for the 2D model. 

Finally, both in the transient and in the steady states 
the currents have small periodic fluctuations determined 
by the permanent transitions of electrons between the 
states in the sample and the states in the leads and 
back.l^They are best seen in Fig. 2(a). Such fluctuations 
have also been obtained very recently by Kurth et at. us- 
ing combination of the non-equilibrium Green's functions 
and the time dependent density- functional theory of the 
Coulomb interaction.!^ 



B. 2D lattice 



case, Eq. (15), except that now the site indices are two- 
dimensional, i.e. i = (ix^iy) and i' = 

The two contact sites are chosen at diagonally opposite 
corners of the sample. The coupling coefficients are cal- 
culated with Eq. (14), like for the ID chain, and depend 



on the wave function of the particular SES at the con- 
tact sites. These coefficients are illustrated in Fig. 5(a). 
The reduced density matrix is calculated using the first 
50 MES. This allowed us to take into account many-body 
configurations with up to 3 electrons. 

In Fig. 5(b) we show the chemical potentials /i^^ for 
the ground and excited states with = 1, 2, and 3 par- 
ticles. At the initial moment to = the system is empty. 
Based on the previous example, the main contribution to 
the currents in the steady states is expected from those 
MES with ground state chemical potentials located inside 
the bias window [fiR, /i^]- One also observes excited con- 
figurations with N particles having chemical potentials 
larger than /i^^x- 
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We show now results for a 2D rectangular lattice with 
12 X 10 sites. For a lattice constant of a = 5 nm this 
sample can be seen as a discrete version of a quantum 
dot of 60 nm x 50 nm. We used the lowest 10 SES of the 
non-interacting sample in the numerical diagonalization 
of the interacting Hamiltonian. This number is sufficient 
to produce convergent results for the first 50 MES for 
an interaction strength U = 0.8. The Coulomb matrix 
elements are calculated in the same way as for the ID 



FIG. 5: (Color online) (a) The coupling amplitudes |Tgn|^ 
for n = 1, .., 5 between single-particle states in the leads with 
momentum q and the lowest 5 single-particle states of the 
isolated dot. (b) The generalized chemical potentials for N - 
particle interacting configurations. The red crosses mark /i^'* 
while the other ones correspond to generalized potentials /x^'* 
related to the i-th excited state of the N particle system. 



In the following we discuss the currents carried by the 
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various many-body states involved in transport. In a first 
series of calculations we selected the chemical potential 
jan = 0.2 and used two values of the chemical potential 
of the left lead /iL = 0.4 and jj^l = 0.6. For juji = 0.2 and 
fiL = 0.4 the bias window contains only the 1-st and the 
2-nd excited configurations with A/" = 1, Fig. 5(b). The 
ground states for N = 1 and N = 2 are instead located 
below and above the bias window, respectively. Conse- 
quently the steady state current is very small. When 
fiL increases to 0.6 the ground state configuration with 
N = 2 enters the bias window and the current increases, 
Fig. 6(a). 

To analyze the transient regime we split the current 
into contributions given by the ground state and excited 
states with 1 electron (see Fig. 6(b)). When ji^ = 0.4 the 
1-st and 2-nd excited state carry currents exceeding the 
current associated to the ground state, which survive all 
the way to the steady state. The current corresponding 
to the 2-nd excited state is smaller than the current of the 
1-st excited state, but comparable to that of the ground 
state. This is explained by the strength of the coupling 
coefficients shown in Fig. 5(a), the 2-nd single-particle 
state being stronger coupled to the leads. The remaining 
higher excited states give oscillating and fast decaying 
transient currents. In Fig. 6(c) /i^ = 0.6 and therefore 
higher excited states enter the bias window; their tran- 
sient currents are still decaying but at a smaller rate. 
Comparing with Fig. 6(a) it in clear that the transient 
regime is dominated by excited states. 

Next we discuss currents associated with states having 
2 and 3 electrons. We keep now fixed jiR = 0.35 and 
again increase fiL starting with 0.6. Fig. 7(a) shows the 
total currents in the left lead for = 2 and = 3. As 
the bias increases the transient currents are enhanced, 
but they become comparable as the system approaches 
the steady state. In Fig. 7(b) the total current on three 
particle states shows a different behavior: the steady 
states value increases drastically when jj^l nioves up. To 
explain this one can look again at the diagram of the 
chemical potentials, Fig. 5(b). At /i^ = 0.6 the 3-particle 
configurations are above the bias window and as such 
they contribute less to the current. In contrast, as /i^ 
increases the ground state configuration with N = 3 en- 
ters the bias window, the window is closer to the excited 
states, and thus the total current increases. Actually, for 
jiiL = 0.8 and 0.9 the current for N = 3 does not reach 
the steady state in the time interval considered. 

Now we look at the contribution of the excited states 
with N = 2 for two cases, /iL = 0.6 and jul = 0.9. Again, 
the inspection of the diagram in Fig. 5(b) predicts the 
results of Fig. 8. When /j^l =0.6 there is just one excited 
configuration within the bias window, in addition to the 
ground state. In Fig. 8(a) we see that in the steady state 
these two configurations give significant contributions to 
the current, whereas the higher excited states play a role 
only in the transient regime. Fig. 8(b) shows that at 
jiiL = 0.9 the currents of the excited states and of the 
ground state are decreasing, some of them reaching even 
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FIG. 6: (a) The total currents in the left and right leads for 
IJiL — 0.6 and iil — 0.4, while keeping iir — 0.2. (b) The 
partial currents in the left lead for single-particle states when 
11 L — 0-4 and iir = 0.2. (c) The partial currents in the left 
lead for single-particle states when fiL = 0.6 and fiR = 0.2 



negative values towards the steady state. This happens 
because the bias window includes now the ground state 
with = 3 and the excited states with N = 2 deplete in 
the favor of the ground state. 

The sign of the current carried by states with N par- 
ticles depends on the placement of the corresponding 
ground state chemical potential relatively to the bias win- 
dow. For example if we fix /i/, = 1.5 and jiR = 0.65 we 
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FIG. 7: (a) The total current in the left lead carried by all 
many-body configurations with N — 2^ for increasing values 
of /iL (i e. 0.6,0.7,0.8 and 0.9) and ^lr = 0.2. (b) The same 
for AT = 3. 



FIG. 8: (a) The total current in the left lead carried by all 
many-body configurations with A/' = 2at/iL=0.6. (b) The 
same for iil — 0.9. Other parameters fiR = 0.35. 



C. Parabolic quantum wire 



obtain fi^ < f^L- Fig. 9(a) shows the A/'-particle cur- 
rents when the sample initially contains two electrons in 
the ground state. This initial state evolves faster to the 
steady state than the empty system. While for = 3 
the current in the left lead is positive, for both N = 2 
and N = 1 the currents are negative. The charge re- 
siding on each A^-particle state and the total charge are 
shown in Fig. 9(b). Since single-particle configurations 
are unlikely their occupation vanishes. The total charge 
accumulated on the N = 3 states increases up to 2, while 
the total charge on the N = 2 states decreases from 2 to 
0.75. The sign of the current for = 2 becomes pos- 
itive when jiji is lowered to 0.2, Fig. 9(c), and exceeds 
the current carried by the states with N = 3. This is 
because the 1-st and the 2-nd SES practically determine 
the ground state with two electrons and thus /i2^\ and 
also because the 1-st SES is strongly coupled to the leads. 
However, the current with A/" = 1 is still negative. 



In this subsection we apply the GME with Coulomb 
interaction to describe the transport through a short 
quantum wire of length = 300 nm with a parabolic 
confinement in the ^/-direction perpendicular to the di- 
rection of transport. The contact ends of the isolated 
wire at ±La^/2 are described by hard walls. This is 
now a continuous model, where a large functional ba- 
sis is used to expand the eigenfunctions of the system 
in. In a similar manner we use a functional basis with 
complete truncated sets of continuous and discrete func- 
tions to expand the eigenfunctions of the semi-infinite 
parabolic leads in. To show that we can describe the com- 
bined geometrical effects imposed on the system by it's 
geometry and an external perpendicular magnetic field 
we place the quantum wire is in an external magnetic 
field of strength 1.0 T. The characteristic confinement 
energy is given by HVLq = 1.0 meV. We assume GaAs 
parameters with m* = 0.067me, = 12.4 meV. The 
magnetic length modified by the parabolic confinement 
is = \/h/(m*rtw)^ with = ^Lq ^ ujI. and the 
cyclotron frequency uOc = eB/{m^c). At B = 1.0 T, 
= 23.87 nm. The semi-infinite leads having the same 



9 



\l^=^ .5, Hr=0.65, N=2 




1 


0.8 






0.6 


rren 


0.4 






O 






0.2 









-0.2 




2 




1.8 




1.6 




1.4 


Q) 


1.2 


D) 


1 


CO 






o 


0.8 




0.6 




0.4 




0.2 








^Il=1 .5, ^ipi=0.65, N=2 



(b) 



N3 
NtOt 



V 



10 20 30 

Time (ps) 

^Il=1.2, ^Fi=0-20. N=0 



40 



50 




20 30 
Time (ps) 

^i|_=1.2, ^iR=0.20, N=0 



50 









N^di 






✓ 




/ 

" ; / 

; / 









10 



20 30 
Time (ps) 



40 



50 



FIG. 9: (a) The total current in the left lead carried by N- 
particle states and the total charge, for /j^l — 1.5 and for 
IJiR = 0.65. (b) The occupation number of the A^-particle 
states, (c) The total current in the left lead carried by N- 
particle states ioi fiL = 1-2 and /i^ = 0.2. (d) The occupation 
number of the A/'-particle states and the total charge. 



parabolic confinement and being subject to the same ex- 
ternal perpendicular magnetic field have a continuous en- 
ergy spectrum with discrete Landau sub-bands. 

The Coulomb potential in Eq. Q in the 2D wire is 
described by 



u{y — y') = 



(16) 



with the small convergence parameter (rj/a^) = 0.01 
to facilitate the two-dimensional numerical integration 
needed for the matrix elements Q. 

After the GME ([9| has been transformed to the in- 
teracting many-electron basis by the unitary transfor- 
mation obtained by the diagonalization of Hs ^ we 
truncate the RDO ^ to 32 MES. For the bias range 
0.0 < A/j. = fiL - /J^R < 1.7 meV used here 10 SES are 
sufficient to obtain these lowest 32 states with good accu- 
racy. We will be omitting singly occupied states of high 
energy that should not be relavant for the parameters 
here. The natural strength of the Coulomb interaction 
will only give us MES that are occupied by one or two 
electrons in the energy range to 6 meV covered by the 
32 MES. 

Since in the partitioning approach [Hs^Hl] = we 
have to construct T^^ as a non-local overlap of (pn and 
ipq'^ on the contact regions C^, I = i?:— 



drdr' (^;(r)^„(r)5^„(r, r') + h.c.) . (17) 



. \E„ 

X exp 



xr-Si{y-y'f 



(18) 



As before Sqi is the energy spectrum of lead /, and En 
is the energy of the SES numbered by n in the quan- 
tum wire. The quantum number q for the states in leads 
represents both the discrete Landau band number and a 
continuous quantum number that can be related to the 
momentum of a particular state. Here we use the pa- 
rameters Sialj 820?^ 0.25, A^^ 0.25 meV, and 
g^R = 40 meV for B = 1.0 T. The domain of the over- 
lap integral for the leads is ±2a^(; into the lead or the 
system for x and x' from each end of the wire at ±La^/2 
and between ±4a^ for y and y' ^ see Ref. (|8]) for an exact 
definition. All the SES will be coupled to the leads, but 
the coupling strength will depend on the character of the 
SES, whether it is an edge- or bulk state and other finer 
geometrical details that is brought about by the magnetic 
field. 

The right chemical potential jiiji is held at 1.4 meV 
and the transport properties are calculated for different 



values of the bias A/i by varying /i^,. Figure [TO] compares 
the total occupation of all one- electron and two-electron 
MES for the interacting system at two different values of 
the bias. At, Aja = 0.2 meV we see that almost solely 
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FIG. 10: (Color online) The total charge residing in one- and 
two-electron states as a function of time for two different val- 
ues of the bias A/i. B 1.0 T, = 300 nm, hQo = 1.0 
meV. 
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FIG. 11: The total steady state current for interacting 10 
SES, and the total charge at t = 270 ps. for different values 
of the bias A/i. B = 1.0 T, = 300 nm, hQo = 1.0 meV. 



32 MES selected here include no three-electron or MES 
with higher number of electrons. It should be mentioned 
here that a different choice of the right bias jiiR can result 
in the system charging faster and thus at the same time 
the total current through it being smaller. This comes 
from the fact that the states have a different coupling to 
the leads and the time range shown here is very much in 
the transient- or it's long exponential decay regime. 



Figure 12 displaying the current in the right lead gives 
an idea how the Coulomb blocking plateau appears after 
the transition regime. The transition regime where the 



one-electron states are occupied, while for A/i = 1.2 meV 
initially it is likely to have one-electron states occupied, 
but very soon the occupation of the two-electron states 
becomes as probable with the likelihood of the occupa- 
tion of the one-electron states fast reducing with time. 
We also have to admit here that even though the steady 
state value of the total current through the system can 
be deduced by the values of the current at 270 ps, the 
charging of the system takes much longer time, since we 
are using here a very weak coupling to the leads that 
mimics a tunneling regime. 

If we now use the average value of the current in the left 
and right leads at t = 270 ps as a measure of the steady 
state current we get the information displayed in Fig. 
[TT} where the steady state value of the current is shown 
for the interacting system as a function of the bias and 
compared to the charge in the system. We have a clear 
Coulomb blocking in the interacting system. In the case 
of a non-interacting system the lack of a gap between the 
one- and two electron MES and a strong mixing of the 
energy regimes of two- and three-electron states the two- 
electron plateau only appears as a small shoulder. The 




FIG. 12: The total current in the right lead for interacting 
and non- interacting 10 SES as a function of the bias A/i and 
time. 5 = 1.0 T, La. = 300 nm, hQo = 1.0 meV. 

right current goes negative, i. e. where it supplies charge 
to the system is partially truncated from the figure. 
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IV. SUMMARY AND CONCLUSIONS 

We calculated time-dependent currents in open meso- 
scopic systems composed by a sample attached to two 
semi-infinite leads, by solving the generalized master 
equation for the reduced density operator acting in the 
Fock space of the sample. This is the natural frame- 
work for including the Coulomb electron-electron inter- 
action in the sample, which is the main achievement of 
this work. The Coulomb interaction is treated in the 
spirit of the exact diagonalization method, i. e. in a pure 
many-body manner. The interacting many-body states 
of the sample are expanded in the basis of non-interacting 
"bit-string" states with unspecified number of electrons. 
We believe our method is a viable alternative to a recent 
appro ach bas ed on a time-dependent density-functional 
model.^ ^^ l ^^ l We used three sample models, a short ID 
wire with 5 sites, but also a larger 2D lattice with 120 
sites and a continuous model, whereas the cited group 
used much smaller samples even with no structure. 

Indeed, due to computational limitations we could use 
only a restricted, effective number of many-body states in 
the GME, between 30-50 depending on the model, from 
the bottom of the energy spectrum. We chose the bias 
window [iir^iil] and the strength of the sample-leads 
coupling parameters Vr^l such that only the effective 
states contribute to the transport of electrons through 
the sample, whereas the states with higher energy are 
unreachable by the electrons. Consequently the number 
of electrons in the sample can be only up to 3 or 4. 

We could calculate the contribution to the charge and 
currents in the sample and in the leads respectively, cor- 
responding to any particular many-body state. We use 
the ID chain as a toy model to emphasize the dominant 
role of the excited states in the transient regime and the 
onset of the Coulomb blockade in the steady state. A 
similar ID model with 4 sites ID has been considered 
recently by Myohanen et al¥^ 

As shown also in our previous works on time-dependent 



transport in non-interacting systems the GME method 
includes information on the energy structure of the sam- 
ple, but also on the geometrical properties reflected in 
the wave functions and sample- lead contacts) ^ * ^ * ^^ * Here 
we illustrate these aspects, in the interacting case, for two 
nanosystems: a two-dimensional quantum dot described 
by a lattice Hamiltonian and a short parabolic quantum 
wire. The time-dependent occupation of speciflc many- 
body states was thoroughly analyzed, for different values 
of the chemical potentials of the leads. It turned out 
that the excited states with N electrons contribute to 
the steady state currents if the ground state conflguration 
with + 1 particles is not available for transport. How- 
ever, if fi^^ < fiR and at the same time /i^^^ lies within 
the bias window the excited states with N particles 
are active only in the transient regime and become de- 
populated in the steady state regime. This behavior is of 
interest in the excited-state spectroscopy experiments. HI 
To our knowledge the time-dependent currents associated 
to excited states have not been discussed theoretically so 
far. 
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